### Set Directory
try(setwd(dirname(rstudioapi::getActiveDocumentContext()$path)))
rm(list=ls())

### Load Package
library(stargazer)
library(readxl)
library(dplyr)
library(tidyr)
library(Hmisc)

### Load Data
data = read.csv('replication dataset.csv')

# Restrict sample to potential targets
# that have economic sanctions experienced
# threat and/or imposition during 1981-2005
data_main = data %>%
  filter(year %in% 1981:2005) %>%
  arrange(ccode,year) %>%
  group_by(ccode) %>%
  mutate(ever = ifelse(mean(threat)>0|mean(imposition)>0,1,0)) %>%
  ungroup() %>%
  filter(ever==1)

### t-test of Dependent Variables for Human Rights vs. Non-human Rights Sanctions

# Figure OA2
data_main = data_main %>%
  mutate(sanctioncat = ifelse(hrdummy==1,1,
                           ifelse(nonhrdummy==1,2,3))) %>%
  mutate(sanctioncat = factor(sanctioncat,labels=c('hr sanctions','nonhr sanctions','none')))

par(mfrow=c(2,2),mar=c(4,3,3,2))
boxplot(v2juncind.rev.lag~sanctioncat,data_main,
        main='A. Lower Court Dependence',xlab='',ylab='')
boxplot(v2x_corr.lag~sanctioncat,data_main,
        main='B. Political Corruption',xlab='',ylab='')
boxplot(v2mecenefm.rev.lag~sanctioncat,data_main,
        main='C. Government Censorship',xlab='',ylab='')
boxplot(lphysintrev~sanctioncat,data_main,
        main='D. Physical Integrity',xlab='',ylab='')

t.test(v2juncind.rev.lag~sanctioncat,data_main[data_main$sanctioncat!="none",])
t.test(v2x_corr.lag~sanctioncat,data_main[data_main$sanctioncat!="none",])
t.test(v2mecenefm.rev.lag~sanctioncat,data_main[data_main$sanctioncat!="none",])
t.test(lphysintrev~sanctioncat,data_main[data_main$sanctioncat!="none",])

### First Stage OLS
m1 = lm(v2juncind.rev~v2juncind.rev.lag+hrdummy+nonhrdummy+polconiii+yr_sch_sec+open.log, data=data_main)
m2 = lm(v2x_corr~v2x_corr.lag+hrdummy+nonhrdummy+polconiii+yr_sch_sec+open.log, data=data_main)
m3 = lm(v2mecenefm.rev~v2mecenefm.rev.lag+hrdummy+nonhrdummy+polconiii+yr_sch_sec+open.log, data=data_main)

stargazer(m1,m2,m3,type='text')

### Second Stage Bootstrap

# Function to get fitted value
hat = function(d,m){
  X = d %>%
    select(names(coef(m))[-1])
  X = cbind(1,X)
  fit = as.matrix(X)%*%coef(m)
  return(fit)
}

# leave out observations with any missing data
d = data_main %>%
  select(v2juncind.rev.lag,v2x_corr.lag,v2mecenefm.rev.lag,
         hrdummy,nonhrdummy,polconiii,yr_sch_sec,open.log,
         tphysintrev,physintrev,
         gdppclog,civilwar,interwar,
         polity2,exrec,exconst,polcomp) %>%
  drop_na()

### Polity IV Component

# EXREC,EXCONST,POLCOMP
boot = function(x){
  d1 = d[sample(1:nrow(d),replace=T),]
  df1 =d1 %>%
    mutate(fit1 = hat(d1,m1)) %>%
    mutate(fit2 = hat(d1,m2)) %>%
    mutate(fit3 = hat(d1,m3))
  
  mb1 = lm(tphysintrev~physintrev+fit1+fit2+fit3
           +gdppclog+exrec+exconst+polcomp+civilwar+interwar,
           data=df1)
  
  beta1 = coef(mb1)
  
  return(beta1)
}

# Bootstrapping (N=1000 iterations)
set.seed(2020)
N=1000
betas = replicate(N,boot())

# 95% Confidence Interval
result = t(apply(betas,1,function(x){
  coef = mean(x)
  lwr = sort(x)[floor(N*.025)]
  upr = sort(x)[floor(N*.975)]
  lwr2 = sort(x)[floor(N*.05)]
  upr2 = sort(x)[floor(N*.95)]
  return(cbind(coef,lwr,upr,lwr2,upr2))
}))

result = as.data.frame(result)

colnames(result) = c('coef','lwr','upr','lwr2','upr2')
result$id = 1:nrow(result)

# Table OA9
round(result[,1:3],3)
nrow(d)

# Coefficient Plot
dev.off()
par(mar=c(5,2,2,2))
plot(id~coef,data=result[3:5,],ylim=c(6,2),xlim=c(min(result[3:5,-6]),max(result[3:5,-6])),yaxt='n',xlab='Coefficient',ylab='',pch=c(15,16,17))
abline(v=0,lty='dashed')
for(i in 3:5){
  lines(y=c(i,i),x=c(result[i,2],result[i,3]))
}
text(y=result[3:5,6]-0.2,x=result[3:5,1],labels=c('Lower Court Dependence','Political Corruption','Government Censorship'),pos=4)

# EXREC
boot = function(x){
  d1 = d[sample(1:nrow(d),replace=T),]
  df1 =d1 %>%
    mutate(fit1 = hat(d1,m1)) %>%
    mutate(fit2 = hat(d1,m2)) %>%
    mutate(fit3 = hat(d1,m3))
  
  mb1 = lm(tphysintrev~physintrev+fit1+fit2+fit3
           +gdppclog+exrec+civilwar+interwar,
           data=df1)
  
  beta1 = coef(mb1)
  
  return(beta1)
}

# Bootstrapping (N=1000 iterations)
set.seed(2020)
N=1000
betas = replicate(N,boot())

# 95% Confidence Interval
result = t(apply(betas,1,function(x){
  coef = mean(x)
  lwr = sort(x)[floor(N*.025)]
  upr = sort(x)[floor(N*.975)]
  lwr2 = sort(x)[floor(N*.05)]
  upr2 = sort(x)[floor(N*.95)]
  return(cbind(coef,lwr,upr,lwr2,upr2))
}))

result = as.data.frame(result)

colnames(result) = c('coef','lwr','upr','lwr2','upr2')
result$id = 1:nrow(result)

# Table OA10
round(result[,1:3],3)
nrow(d)

# Coefficient Plot
par(mar=c(5,2,2,2))
plot(id~coef,data=result[3:5,],ylim=c(6,2),xlim=c(min(result[3:5,-6]),max(result[3:5,-6])),yaxt='n',xlab='Coefficient',ylab='',pch=c(15,16,17))
abline(v=0,lty='dashed')
for(i in 3:5){
  lines(y=c(i,i),x=c(result[i,2],result[i,3]))
}
text(y=result[3:5,6]-0.2,x=result[3:5,1],labels=c('Lower Court Dependence','Political Corruption','Government Censorship'),pos=4)

# EXCONST
boot = function(x){
  d1 = d[sample(1:nrow(d),replace=T),]
  df1 =d1 %>%
    mutate(fit1 = hat(d1,m1)) %>%
    mutate(fit2 = hat(d1,m2)) %>%
    mutate(fit3 = hat(d1,m3))
  
  mb1 = lm(tphysintrev~physintrev+fit1+fit2+fit3
           +gdppclog+exconst+civilwar+interwar,
           data=df1)
  
  beta1 = coef(mb1)
  
  return(beta1)
}

# Bootstrapping (N=1000 iterations)
set.seed(2020)
N=1000
betas = replicate(N,boot())

# 95% Confidence Interval
result = t(apply(betas,1,function(x){
  coef = mean(x)
  lwr = sort(x)[floor(N*.025)]
  upr = sort(x)[floor(N*.975)]
  lwr2 = sort(x)[floor(N*.05)]
  upr2 = sort(x)[floor(N*.95)]
  return(cbind(coef,lwr,upr,lwr2,upr2))
}))

result = as.data.frame(result)

colnames(result) = c('coef','lwr','upr','lwr2','upr2')
result$id = 1:nrow(result)

# Table OA11
round(result[,1:3],3)
nrow(d)

# Coefficient Plot
par(mar=c(5,2,2,2))
plot(id~coef,data=result[3:5,],ylim=c(6,2),xlim=c(min(result[3:5,-6]),max(result[3:5,-6])),yaxt='n',xlab='Coefficient',ylab='',pch=c(15,16,17))
abline(v=0,lty='dashed')
for(i in 3:5){
  lines(y=c(i,i),x=c(result[i,2],result[i,3]))
}
text(y=result[3:5,6]-0.2,x=result[3:5,1],labels=c('Lower Court Dependence','Political Corruption','Government Censorship'),pos=4)

# POLCOMP
boot = function(x){
  d1 = d[sample(1:nrow(d),replace=T),]
  df1 =d1 %>%
    mutate(fit1 = hat(d1,m1)) %>%
    mutate(fit2 = hat(d1,m2)) %>%
    mutate(fit3 = hat(d1,m3))
  
  mb1 = lm(tphysintrev~physintrev+fit1+fit2+fit3
           +gdppclog+polcomp+civilwar+interwar,
           data=df1)
  
  beta1 = coef(mb1)
  
  return(beta1)
}

# Bootstrapping (N=1000 iterations)
set.seed(2020)
N=1000
betas = replicate(N,boot())

# 95% Confidence Interval
result = t(apply(betas,1,function(x){
  coef = mean(x)
  lwr = sort(x)[floor(N*.025)]
  upr = sort(x)[floor(N*.975)]
  lwr2 = sort(x)[floor(N*.05)]
  upr2 = sort(x)[floor(N*.95)]
  return(cbind(coef,lwr,upr,lwr2,upr2))
}))

result = as.data.frame(result)

colnames(result) = c('coef','lwr','upr','lwr2','upr2')
result$id = 1:nrow(result)

# Table OA 12
round(result[,1:3],3)
nrow(d)

# Coefficient Plot
par(mar=c(5,2,2,2))
plot(id~coef,data=result[3:5,],ylim=c(6,2),xlim=c(min(result[3:5,-6]),max(result[3:5,-6])),yaxt='n',xlab='Coefficient',ylab='',pch=c(15,16,17))
abline(v=0,lty='dashed')
for(i in 3:5){
  lines(y=c(i,i),x=c(result[i,2],result[i,3]))
}
text(y=result[3:5,6]-0.2,x=result[3:5,1],labels=c('Lower Court Dependence','Political Corruption','Government Censorship'),pos=4)

### Correlation among Repression and Polity Variables

# Table OA13
rcorr(as.matrix(data_main[,c('physintrev','v2juncind.rev','v2x_corr','v2mecenefm.rev','polity2','exrec','exconst','polcomp')]))

### Summary Statistics

stargazer(as.data.frame(data_main %>%
                          select(physint,
                                 v2juncind.rev,
                                 v2x_corr,
                                 v2mecenefm.rev,
                                 hrdummy,
                                 nonhrdummy,
                                 polconiii,
                                 yr_sch_sec,
                                 open.log,
                                 gdppclog,
                                 polity2,
                                 exrec,
                                 exconst,
                                 polcomp,
                                 civilwar,
                                 interwar)),type='html',out='summary_pol.doc')
